require(MCMCpack)

#### prepare data ####
elite.2012 <- read.csv("2012UTASP20150910.csv")
# missing values
elite.2012[elite.2012 == 66 | elite.2012 == 99] <- NA
# leave MPs only
winner.2012 <- subset(elite.2012, RESULT > 0)
# data frame including issue items
issue.elite <- cbind(winner.2012[, 38:72])
# labels of issue items
issue.names <- c("Constitutional revision", "Increasing defensive power", 
                 "Pre-emptive defense", "Permanent member of the UNSC", 
                 "Pressure on North Korea", "Collective self-defense", 
                 "Small government", "Public works", 
                 "Increasing public spending", "Increasing consumption tax", 
                 "Consumption tax over 10%", "No pension for the rich", 
                 "Trans-Pacific Partnership", "Inflation target", 
                 "Security before privacy", "Foreign suffrage", 
                 "Immigrant workers", "Moral education", 
                 "Reactivating nuclear plants", 
                 "Voluntarily accepting debris", 
                 "Strengthening HR", "Prime minister election", 
                 "The doshu system", "Local consumption tax", 
                 "SNTV", "Decreasing Diet members", 
                 "US before Asia", "Improving competitive power", 
                 "Protectionism", "Environmental protection", 
                 "Traditional family", "Zero nuclear energy by 2030", 
                 "Free political donations", "Arena legislature", 
                 "Party discipline")

#### estimate the graded IRT model ####
## one-dimensional model
# sign constraints
one.d.constraints <- list(Q5_2 = list(2, "-"), Q5_6 = list(2, "-"))
# run three chains
one.d.result.1 <- MCMCordfactanal(~., factors = 1, 
                                 lambda.constraints = one.d.constraints, 
                                 data = issue.elite, 
                                 burnin = 5000, mcmc = 100000, thin = 25, L0 = 0.01, 
                                 tune = 0.2, seed = 101, store.scores = TRUE)
one.d.result.2 <- MCMCordfactanal(~., factors = 1, 
                                 lambda.constraints = one.d.constraints, 
                                 data = issue.elite, 
                                 burnin = 5000, mcmc = 100000, thin = 25, L0 = 0.01, 
                                 tune = 0.2, seed = 102, store.scores = TRUE)
one.d.result.3 <- MCMCordfactanal(~., factors = 1, 
                                 lambda.constraints = one.d.constraints, 
                                 data = issue.elite, 
                                 burnin = 5000, mcmc = 100000, thin = 25, L0 = 0.01, 
                                 tune = 0.2, seed = 103, store.scores = TRUE)
# convergence check
one.d.result <- mcmc.list(one.d.result.1, one.d.result.2, one.d.result.3)
which(gelman.diag(one.d.result, multivariate = FALSE)[[1]][, 1] > 1.1)
# combine three chains
one.d.result.matrix <- rbind(as.matrix(one.d.result.1), as.matrix(one.d.result.2), 
                             as.matrix(one.d.result.3))
# save MCMC samples
# save(one.d.result.matrix, file = "one_d_result_matrix.Rdata")
# load MCMC samples
# load("one_d_result_matrix.Rdata")
# discrimination parameters
one.d.loadings <- array(NA, c(12000, 35, 2))
one.d.loadings[, , 1] <- one.d.result.matrix[, seq(1, 69, 2)]
one.d.loadings[, , 2] <- one.d.result.matrix[, seq(2, 70, 2)]
# threshold parameters
one.d.thresholds <- array(NA, c(12000, 35, 3))
one.d.thresholds[, , 1] <- one.d.result.matrix[, 71:105]
one.d.thresholds[, , 2] <- one.d.result.matrix[, 106:140]
one.d.thresholds[, , 3] <- one.d.result.matrix[, 141:175]
# latent trait parameters
one.d.scores <- array(NA, c(12000, nrow(issue.elite), 2))
one.d.scores[, , 1] <- 1
one.d.scores[, , 2] <- one.d.result.matrix[, 176:655]

## two-dimensional model
# sign constraints
two.d.constraints <- list(Q5_2 = list(2, "-"), Q5_2 = list(3, 0), 
                          Q5_6 = list(2, "-"), Q5_6 = list(3, 0), 
                          Q5_22 = list(2, 0), Q5_22 = list(3, "+"), 
                          Q5_24 = list(2, 0), Q5_24 = list(3, "+"))
# run three chains
two.d.result.1 <- MCMCordfactanal(~., factors = 2, 
                                 lambda.constraints = two.d.constraints, 
                                 data = issue.elite, 
                                 burnin = 5000, mcmc = 100000, thin = 25, L0 = 0.01, 
                                 tune = 0.2, seed = 201, store.scores = TRUE)
two.d.result.2 <- MCMCordfactanal(~., factors = 2, 
                                 lambda.constraints = two.d.constraints, 
                                 data = issue.elite, 
                                 burnin = 5000, mcmc = 100000, thin = 25, L0 = 0.01, 
                                 tune = 0.2, seed = 202, store.scores = TRUE)
two.d.result.3 <- MCMCordfactanal(~., factors = 2, 
                                 lambda.constraints = two.d.constraints, 
                                 data = issue.elite, 
                                 burnin = 5000, mcmc = 100000, thin = 25, L0 = 0.01, 
                                 tune = 0.2, seed = 203, store.scores = TRUE)
# convergence check
two.d.result <- mcmc.list(two.d.result.1, two.d.result.2, two.d.result.3)
which(gelman.diag(two.d.result, multivariate = FALSE)[[1]][, 1] > 1.1)
# combine three chains
two.d.result.matrix <- rbind(as.matrix(two.d.result.1), as.matrix(two.d.result.2), 
                             as.matrix(two.d.result.3))
# save MCMC samples
# save(two.d.result.matrix, file = "two_d_result_matrix.Rdata")
# load MCMC samples
# load("two_d_result_matrix.Rdata")
# difficulty and discrimination parameters
two.d.loadings <- array(NA, c(12000, 35, 3))
two.d.loadings[, , 1] <- two.d.result.matrix[, c(1, 4, 6, 9, 12, 15, seq(17, 62, 3), 64, 67, seq(69, 99, 3))]
two.d.loadings[, , 2] <- cbind(two.d.result.matrix[, c(2, 5, 7, 10, 13, 16, seq(18, 60, 3))], 0, 
                               two.d.result.matrix[, 65], 0, two.d.result.matrix[, seq(70, 100, 3)])
two.d.loadings[, , 3] <- cbind(two.d.result.matrix[, 3], 0, two.d.result.matrix[, c(8, 11, 14)], 0, 
                               two.d.result.matrix[, c(seq(19, 61, 3), 63, 66, seq(68, 101, 3))])
# threshold parameters
two.d.thresholds <- array(NA, c(12000, 35, 3))
two.d.thresholds[, , 1] <- two.d.result.matrix[, 102:136]
two.d.thresholds[, , 2] <- two.d.result.matrix[, 137:171]
two.d.thresholds[, , 3] <- two.d.result.matrix[, 172:206]
# latent trait parameters
two.d.scores <- array(NA, c(12000, nrow(issue.elite), 3))
two.d.scores[, , 1] <- 1
two.d.scores[, , 2] <- two.d.result.matrix[, seq(207, 1165, 2)]
two.d.scores[, , 3] <- two.d.result.matrix[, seq(208, 1166, 2)]

## three-dimensional model
# sign constraints
three.d.constraints <- list(Q5_2 = list(2, "-"), Q5_2 = list(3, 0), Q5_2 = list(4, 0), 
                            Q5_6 = list(2, "-"), Q5_6 = list(3, 0), Q5_6 = list(4, 0), 
                            Q5_22 = list(2, 0), Q5_22 = list(3, "+"), Q5_22 = list(4, 0), 
                            Q5_24 = list(2, 0), Q5_24 = list(3, "+"), Q5_24 = list(4, 0), 
                            Q5_10 = list(2, 0), Q5_10 = list(3, 0), Q5_10 = list(4, "+"), 
                            Q5_11 = list(2, 0), Q5_11 = list(3, 0), Q5_11 = list(4, "+"))
# run three chains
three.d.result.1 <- MCMCordfactanal(~., factors = 3, 
                                   lambda.constraints = three.d.constraints, 
                                   data = issue.elite, 
                                   burnin = 5000, mcmc = 100000, thin = 25, L0 = 0.01, 
                                   tune = 0.2, seed = 301, store.scores = TRUE)
three.d.result.2 <- MCMCordfactanal(~., factors = 3, 
                                   lambda.constraints = three.d.constraints, 
                                   data = issue.elite, 
                                   burnin = 5000, mcmc = 100000, thin = 25, L0 = 0.01, 
                                   tune = 0.2, seed = 302, store.scores = TRUE)
three.d.result.3 <- MCMCordfactanal(~., factors = 3, 
                                   lambda.constraints = three.d.constraints, 
                                   data = issue.elite, 
                                   burnin = 5000, mcmc = 100000, thin = 25, L0 = 0.01, 
                                   tune = 0.2, seed = 303, store.scores = TRUE)
# convergence check
three.d.result <- mcmc.list(three.d.result.1, three.d.result.2, three.d.result.3)
which(gelman.diag(three.d.result, multivariate = FALSE)[[1]][, 1] > 1.1)
# combine three chains
three.d.result.matrix <- rbind(as.matrix(three.d.result.1), as.matrix(three.d.result.2), 
                               as.matrix(three.d.result.3))
# save MCMC samples
# save(three.d.result.matrix, file = "three_d_result_matrix.Rdata")
# load MCMC samples
# load("three_d_result_matrix.Rdata")
# discrimination parameters
three.d.loadings <- array(NA, c(12000, 35, 4))
three.d.loadings[, , 1] <- three.d.result.matrix[, c(1, 5, 7, 11, 15, 19, seq(21, 33, 4), 35, 
                                                     seq(37, 77, 4), 79, 83, seq(85, 125, 4))]
three.d.loadings[, , 2] <- cbind(three.d.result.matrix[, c(2, 6, 8, 12, 16, 20, 22, 26, 30)], 0, 0, 
                                 three.d.result.matrix[, seq(38, 74, 4)], 0, three.d.result.matrix[, 80], 
                                 0, three.d.result.matrix[, seq(86, 126, 4)])
three.d.loadings[, , 3] <- cbind(three.d.result.matrix[, 3], 0, three.d.result.matrix[, c(9, 13, 17)], 0, 
                                 three.d.result.matrix[, c(23, 27, 31)], 0, 0, 
                                 three.d.result.matrix[, c(seq(39, 75, 4), 78, 81, 84, seq(87, 127, 4))])
three.d.loadings[, , 4] <- cbind(three.d.result.matrix[, 4], 0, three.d.result.matrix[, c(10, 14, 18)], 0, 
                                 three.d.result.matrix[, c(24, 28, 32, 34, seq(36, 76, 4))], 0, 
                                 three.d.result.matrix[, 82], 0, three.d.result.matrix[, seq(88, 128, 4)])
# threshold parameters
three.d.thresholds <- array(NA, c(12000, 35, 3))
three.d.thresholds[, , 1] <- three.d.result.matrix[, 129:163]
three.d.thresholds[, , 2] <- three.d.result.matrix[, 164:198]
three.d.thresholds[, , 3] <- three.d.result.matrix[, 199:233]
# latent trait parameters
three.d.scores <- array(NA, c(12000, nrow(issue.elite), 4))
three.d.scores[, , 1] <- 1
three.d.scores[, , 2] <- three.d.result.matrix[, seq(234, 1671, 3)]
three.d.scores[, , 3] <- three.d.result.matrix[, seq(235, 1672, 3)]
three.d.scores[, , 4] <- three.d.result.matrix[, seq(236, 1673, 3)]

## four-dimensional model
# sign constraints
four.d.constraints <- list(Q5_2 = list(2, "-"), Q5_2 = list(3, 0), Q5_2 = list(4, 0), Q5_2 = list(5, 0), 
                           Q5_6 = list(2, "-"), Q5_6 = list(3, 0), Q5_6 = list(4, 0), Q5_6 = list(5, 0), 
                           Q5_22 = list(2, 0), Q5_22 = list(3, "+"), Q5_22 = list(4, 0), Q5_22 = list(5, 0), 
                           Q5_24 = list(2, 0), Q5_24 = list(3, "+"), Q5_24 = list(4, 0), Q5_24 = list(5, 0), 
                           Q5_10 = list(2, 0), Q5_10 = list(3, 0), Q5_10 = list(4, "+"), Q5_10 = list(5, 0), 
                           Q5_11 = list(2, 0), Q5_11 = list(3, 0), Q5_11 = list(4, "+"), Q5_11 = list(5, 0), 
                           Q5_8 = list(2, 0), Q5_8 = list(3, 0), Q5_8 = list(4, 0), Q5_8 = list(5, "+"), 
                           Q5_9 = list(2, 0), Q5_9 = list(3, 0), Q5_9 = list(4, 0), Q5_9 = list(5, "+"))
# run three chains
four.d.result.1 <- MCMCordfactanal(~., factors = 4, 
                                  lambda.constraints = four.d.constraints, 
                                  data = issue.elite, 
                                  burnin = 5000, mcmc = 100000, thin = 25, L0 = 0.01, 
                                  tune = 0.2, seed = 401, store.scores = TRUE)
four.d.result.2 <- MCMCordfactanal(~., factors = 4, 
                                  lambda.constraints = four.d.constraints, 
                                  data = issue.elite, 
                                  burnin = 5000, mcmc = 100000, thin = 25, L0 = 0.01, 
                                  tune = 0.2, seed = 402, store.scores = TRUE)
four.d.result.3 <- MCMCordfactanal(~., factors = 4, 
                                  lambda.constraints = four.d.constraints, 
                                  data = issue.elite, 
                                  burnin = 5000, mcmc = 100000, thin = 25, L0 = 0.01, 
                                  tune = 0.2, seed = 403, store.scores = TRUE)
# convergence check
four.d.result <- mcmc.list(four.d.result.1, four.d.result.2, four.d.result.3)
which(gelman.diag(four.d.result, multivariate = FALSE)[[1]][, 1] > 1.1)
# combine three chains
four.d.result.matrix <- rbind(as.matrix(four.d.result.1), as.matrix(four.d.result.2), 
                              as.matrix(four.d.result.3))
# save MCMC samples
# save(four.d.result.matrix, file = "four_d_result_matrix.Rdata")
# load MCMC samples
# load("four_d_result_matrix.Rdata")
# discrimination parameters
four.d.loadings <- array(NA, c(12000, 35, 5))
four.d.loadings[, , 1] <- four.d.result.matrix[, c(1, 6, 8, 13, 18, 23, 25, 30, 32, 34, 36, seq(38, 88, 5), 
                                                   90, 95, seq(97, 147, 5))]
four.d.loadings[, , 2] <- cbind(four.d.result.matrix[, c(2, 7, 9, 14, 19, 24, 26)], 0, 0, 0, 0, 
                                four.d.result.matrix[, seq(39, 84, 5)], 0, four.d.result.matrix[, 91], 
                                0, four.d.result.matrix[, seq(98, 148, 5)])
four.d.loadings[, , 3] <- cbind(four.d.result.matrix[, 3], 0, four.d.result.matrix[, c(10, 15, 20)], 0, 
                                four.d.result.matrix[, 27], 0, 0, 0, 0, 
                                four.d.result.matrix[, c(seq(40, 85, 5), 89, 92, 96, seq(99, 149, 5))])
four.d.loadings[, , 4] <- cbind(four.d.result.matrix[, 4], 0, four.d.result.matrix[, c(11, 16, 21)], 0, 
                                four.d.result.matrix[, 28], 0, 0, 
                                four.d.result.matrix[, c(35, 37, seq(41, 86, 5))], 0, four.d.result.matrix[, 93], 
                                0, four.d.result.matrix[, seq(100, 150, 5)])
four.d.loadings[, , 5] <- cbind(four.d.result.matrix[, 5], 0, four.d.result.matrix[, c(12, 17, 22)], 0, 
                                four.d.result.matrix[, c(29, 31, 33)], 0, 0, 
                                four.d.result.matrix[, seq(42, 87, 5)], 0, four.d.result.matrix[, 94], 
                                0, four.d.result.matrix[, seq(101, 151, 5)])
# threshold parameters
four.d.thresholds <- array(NA, c(12000, 35, 5))
four.d.thresholds[, , 1] <- four.d.result.matrix[, 152:186]
four.d.thresholds[, , 2] <- four.d.result.matrix[, 187:221]
four.d.thresholds[, , 3] <- four.d.result.matrix[, 222:256]
# latent trait parameters
four.d.scores <- array(NA, c(12000, nrow(issue.elite), 5))
four.d.scores[, , 1] <- 1
four.d.scores[, , 2] <- four.d.result.matrix[, seq(257, 2173, 4)]
four.d.scores[, , 3] <- four.d.result.matrix[, seq(258, 2174, 4)]
four.d.scores[, , 4] <- four.d.result.matrix[, seq(259, 2175, 4)]
four.d.scores[, , 5] <- four.d.result.matrix[, seq(260, 2176, 4)]

## model comparison
# function to compute WAIC
information.criterion <- function(y, loading, threshold, score) {
  likelihood <- matrix(NA, nrow(loading), ncol(y) * nrow(y))
  for (i in 1:nrow(y)) {
    for (j in 1:ncol(y)) {
      if (is.na(y[i, j])) {
        likelihood[, (ncol(y) * (i - 1) + j)] <- NA
      } else if (y[i, j] == 1) {
        likelihood[, (ncol(y) * (i - 1) + j)] <- pnorm(-rowSums(loading[, j, ] * score[, i, ]))
      } else if (y[i, j] == 2) {
        likelihood[, (ncol(y) * (i - 1) + j)] <- pnorm(threshold[, j, 1] - rowSums(loading[, j, ] * score[, i, ])) - 
          pnorm(-rowSums(loading[, j, ] * score[, i, ]))
      } else if (y[i, j] == 3) {
        likelihood[, (ncol(y) * (i - 1) + j)] <- pnorm(threshold[, j, 2] - rowSums(loading[, j, ] * score[, i, ])) - 
          pnorm(threshold[, j, 1] - rowSums(loading[, j, ] * score[, i, ]))
      } else if (y[i, j] == 4) {
        likelihood[, (ncol(y) * (i - 1) + j)] <- pnorm(threshold[, j, 3] - rowSums(loading[, j, ] * score[, i, ])) - 
          pnorm(threshold[, j, 2] - rowSums(loading[, j, ] * score[, i, ]))
      } else if (y[i, j] == 5) {
        likelihood[, (ncol(y) * (i - 1) + j)] <- 1 - pnorm(threshold[, j, 3] - rowSums(loading[, j, ] * score[, i, ]))
      } 
    }
  }
  WAIC <- -2 * (sum(log(colMeans(likelihood, na.rm = TRUE)), na.rm = TRUE) - 
                  sum(apply(log(likelihood), 2, var, na.rm = TRUE), na.rm = TRUE))
  WAIC
}

# compute WIACs (Appendix, page 7, lines 8-10 and footnote 2)
round(information.criterion(issue.elite, one.d.loadings, one.d.thresholds, one.d.scores))
round(information.criterion(issue.elite, two.d.loadings, two.d.thresholds, two.d.scores))
round(information.criterion(issue.elite, three.d.loadings, three.d.thresholds, three.d.scores))
round(information.criterion(issue.elite, four.d.loadings, four.d.thresholds, four.d.scores))

# issue items related to the third and fourth dimensions (Appendix, footnote 2)
two.d.loadings.table <- apply(two.d.loadings, c(2, 3), mean)
three.d.loadings.table <- apply(three.d.loadings, c(2, 3), mean)
four.d.loadings.table <- apply(four.d.loadings, c(2, 3), mean)
rownames(two.d.loadings.table) <- issue.names
rownames(three.d.loadings.table) <- issue.names
rownames(four.d.loadings.table) <- issue.names
colnames(two.d.loadings.table) <- c("difficulty", "discrim.1", "discrim.2")
colnames(three.d.loadings.table) <- c("difficulty", "discrim.1", "discrim.2", "discrim.3")
colnames(four.d.loadings.table) <- c("difficulty", "discrim.1", "discrim.2", "discrim.3", "discrim.4")
round(two.d.loadings.table, 2)
round(three.d.loadings.table, 2)
round(four.d.loadings.table, 2)

#### estimation results of discrimination parameters (Figure 1) ####
## compute posterior means and 95% HPD intervals
# posterior means
two.d.discrim.1.mean <- colMeans(two.d.loadings[, , 2])
two.d.discrim.2.mean <- colMeans(two.d.loadings[, , 3])
# 95% HPD intervals
two.d.discrim.1.HPD <- HPDinterval(mcmc(two.d.loadings[, , 2]))
two.d.discrim.2.HPD <- HPDinterval(mcmc(two.d.loadings[, , 3]))
# reverse negative items
two.d.discrim.1.results <- cbind(two.d.discrim.1.mean, two.d.discrim.1.HPD, 0)
two.d.discrim.2.results <- cbind(two.d.discrim.2.mean, two.d.discrim.2.HPD, 0)
for (i in 1:35) {
  if (two.d.discrim.1.results[i, 1] < 0) {
    two.d.discrim.1.results[i, ] <- -two.d.discrim.1.results[i, ]
    two.d.discrim.1.results[i, 4] <- 1
  }
  if (two.d.discrim.2.results[i, 1] < 0) {
    two.d.discrim.2.results[i, ] <- -two.d.discrim.2.results[i, ]
    two.d.discrim.2.results[i, 4] <- 1
  }
}

## draw the figure
rownames(two.d.discrim.1.results) <- rownames(two.d.discrim.2.results) <- issue.names
# descending order
order1 <- order(two.d.discrim.1.results[, 1])
order2 <- order(two.d.discrim.2.results[, 1])
two.d.discrim.1.results <- two.d.discrim.1.results[order1, ]
two.d.discrim.2.results <- two.d.discrim.2.results[order2, ]

tiff("Fig1.tif", width = 6, height = 6, units = "in", pointsize = 6, res = 300)
layout(matrix(1:2, 1, 2))
dotchart(two.d.discrim.1.results[, 1], pch = 20, xlim = c(-0.3, 3.8), 
         main = "Liberal\u2014Conservative", xlab = "Discrimination parameters")
segments(two.d.discrim.1.results[, 2], 1:35, two.d.discrim.1.results[, 3], 1:35)
points(rep(-0.3, 35), 1:35, pch = ifelse(two.d.discrim.1.results[, 4] == 1, 17, NA), cex = 0.6)
abline(v = 0, lty = 3)
dotchart(two.d.discrim.2.results[, 1], pch = 20, xlim = c(-0.3, 1.6), 
         main = "Radical-Reform\u2014Maintenance", xlab = "Discrimination parameters")
segments(two.d.discrim.2.results[ ,2], 1:35, two.d.discrim.2.results[, 3], 1:35)
points(rep(-0.3, 35), 1:35, pch = ifelse(two.d.discrim.2.results[, 4] == 1, 17, NA), cex = 0.6)
abline(v = 0, lty = 3)
dev.off()

#### estimation results of latent traits (Figure 2) ####
## discard MPs who did not answer two-thirds or more of the items
n.NA <- function(x) sum(is.na(x))  # function to compute the sum of non-responses
elite.NA <- ifelse(apply(issue.elite, 1, n.NA) > 11, 1, 0)
elite.score <- list()
elite.score[[1]] <- two.d.scores[, , 2][, elite.NA == 0]
elite.score[[2]] <- two.d.scores[, , 3][, elite.NA == 0]

## party affiliation
elite.party <- winner.2012$PARTY[elite.NA == 0]
party.score <- list()
for (i in 1:7) {
  party.score[[i]] <- list()
  for (j in 1:2) {
    party.score[[i]][[j]] <- elite.score[[j]][, elite.party == i]
  }
}

## draw the figure
colors <- c("#ff2800", "#0041ff", "#35a16b", "#7f878f")
tiff("Fig2.tif", width = 6, height = 6, units = "in", pointsize = 9, res = 300)
par(mar = c(4, 4, 2, 2))
plot(0, 0, type = "n", xlim = c(-4, 2.5), ylim = c(-4, 2.5), 
     xlab = "Liberal\u2014Conservative", ylab = "Radical-Reform\u2014Maintenance")
points(colMeans(party.score[[4]][[1]]), colMeans(party.score[[4]][[2]]), 
       pch = "G", col = colors[4], cex = 0.7)
points(colMeans(party.score[[6]][[1]]), colMeans(party.score[[6]][[2]]), 
       pch = "C", col = colors[4], cex = 0.7)
points(colMeans(party.score[[7]][[1]]), colMeans(party.score[[7]][[2]]), 
       pch = "Y", col = colors[4], cex = 0.7)
points(colMeans(party.score[[2]][[1]]), colMeans(party.score[[2]][[2]]), 
       pch = "L", col = colors[1], cex = 0.7)
points(colMeans(party.score[[1]][[1]]), colMeans(party.score[[1]][[2]]), 
       pch = "D", col = colors[2], cex = 0.7)
points(colMeans(party.score[[5]][[1]]), colMeans(party.score[[5]][[2]]), 
       pch = "R", col = colors[3], cex = 0.7)
text(mean(party.score[[2]][[1]]), mean(party.score[[2]][[2]]), labels = "L", cex = 1.5, font = 2)
text(mean(party.score[[1]][[1]]), mean(party.score[[1]][[2]]), labels = "D", cex = 1.5, font = 2)
text(mean(party.score[[5]][[1]]), mean(party.score[[5]][[2]]), labels = "R", cex = 1.5, font = 2)
legend("bottomright", legend = c("Liberal Democratic Party", "Democratic Party of Japan", 
                                 "Japan Restoration Party", "Clean Government Party", 
                                 "Your Party", "Japanese Communist Party"), 
       pch = c("L", "D", "R", "G", "Y", "C"), col = colors[c(1, 2, 3, 4, 4, 4)], cex = 0.9, pt.cex = 0.7)
dev.off()